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Abstract.  Numerical  methods  for  grids  with  irregular  cells  require  discrete  shape  functions  to  approximate  the  distribution  of 
quantities  across  cells.  For  control- volume  mixed  finite-element  (CVMFE)  methods,  vector  shape  functions  approximate  velocities 
and  vector  test  functions  enforce  a  discrete  form  of  Darcy’s  law.  In  this  paper,  a  new  vector  shape  function  is  developed  for  use  with 
irregular,  hexahedral  cells  (trilinear  images  of  cubes).  It  interpolates  velocities  and  fluxes  quadratically,  because  as  shown  here, 
the  usual  Piola-transformed  shape  functions,  which  interpolate  linearly,  cannot  match  uniform  flow  on  general  hexahedral  cells. 
Truncation-error  estimates  for  the  shape  function  are  demonstrated.  CVMFE  simulations  of  uniform  and  non-uniform  flow  with 
irregular  meshes  show  first-  and  second-order  convergence  of  fluxes  in  the  Lr  norm  in  the  presence  and  absence  of  singularities, 
respectively. 


1.  INTRODUCTION 


For  simulation  of  two-dimensional  (2-D)  flow  in  heterogeneous  porous  media,  it  has  been  shown  that  mixed 
methods,  and  in  particular  the  control-volume  mixed  finite-element  (CVMFE)  method  [1],  can  be  efficient, 
accurate  schemes  to  obtain  the  velocity  field  [1-3].  As  part  of  an  extension  of  CVMFE  to  3-D,  here  we 
develop  and  test  a  3-D  CVMFE  velocity  shape  function  for  irregular  hexahedral  cells,  based  on  covariant 
vectors  for  a  mapping  to  a  unit  cube.  Its  performance  is  evaluated  with  the  L2  norm  of  the  flux  error  in  a 
number  of  test  cases. 

In  3-D  CVMFE,  the  domain  is  discretized  into  potentially  irregular  hexahedral  cells,  allowing  for  the 
modeling  of  complex  hydrogeological  systems.  Vector  shape  functions,  described  subsequently,  interpolate 
the  velocity  over  the  interiors  of  cells,  and  vector  test  functions  are  weighting  factors  in  integrating  the  Darcy 
relation  over  control  volumes  associated  with  cell  faces;  this  can  be  viewed  as  an  error  minimization  in  the 
control- volume  technique  [4].  The  shape  and  test  functions  result  in  discrete  equations  that  can  be  solved  for 
fluxes  at  cell  faces  and  pressures  at  cell  centers.  Shape  functions  interpolate  cell-surface  fluxes,  approximat¬ 
ing  the  velocity  in  the  cell  interior;  a  poor  approximation  degrades  the  solution.  We  propose  a  new  velocity 
shape  function  for  3-D  logically  rectangular  meshes  which,  in  most  cases,  should  provide  a  reasonable  cell- 
velocity  estimate.  This  function  is  based  on  a  second-order  approximation  of  flux  conservation  across  the 
cell,  which  can  be  a  non-linear  interpolation  of  the  flux  and  the  velocity.  For  3-D  flow  on  irregular  meshes, 
we  believe  that  this  function  will  offer  advantages  over  3-D  versions  of  the  2-D  linear  shape  function  in  [1]. 

Sections  2-5  provide  background  on  the  problem  and  the  solution  technique.  In  particular,  section  3 
describes  the  roles  of  the  shape  functions  and  test  functions  in  CVMFE.  Section  4  shows  that  the  flux  of 
a  uniform  flow  field  across  an  arbitrary  cell  is  not  necessarily  linear;  this  affects  the  construction  of  the 
shape  functions.  Section  5  discusses  the  Piola  transformation,  which  determines  standard  shape  functions 
with  lineal-  flux  for  irregular  cells.  Velocity  interpolation  functions,  a  precursor  of  the  shape  functions,  are 
presented  in  section  6,  followed  by  the  shape  functions  themselves  in  section  7.  Sections  8  and  9  present 
results  of  the  velocity  shape  function  for  uniform  and  non-uniform  flow,  respectively.  A  brief  discussion  of 
the  results  concludes  the  paper. 
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Figure  1.  Left:  reference  cube  Q  with  edges  of  unit  length.  Right:  arbitrary  cell  Q  from  discretization  with  vertex  locations  vooo. 
Vioo,  Voio,  voot,  V110,  Viol,  von  and  vm  indicated. 


2.  BASIC  EQUATIONS 

This  paper  applies  CVMFE  to  a  steady  flow  equation  in  a  3-D  domain  Q.\ 

V-q  =  w{x,y,z),  {x,y,z)  E  £1  (1) 

Here  q  is  the  specific  discharge  vector  and  W  is  a  source  term.  Boundary  conditions  on  the  surface  dQ  are 
fluxes  over  c)flj  and/or  specified  pressures  over  dQp.  We  assume  the  Darcy  relation: 

q  =  -K  Vp/fj,  (2) 

where  p  is  the  pressure,  K (x,y,z)  is  the  hydraulic  conductivity  tensor  and  p  is  the  viscosity  (the  gravitational 
potential  is  neglected  from  (2)  for  notational  convenience).  Like  most  mixed  methods,  CVMFE  inverts  the 
hydraulic  conductivity  tensor  in  (2) ,  so  that 

Vp  =  -p  K_1q.  (3) 

The  mixed- method  development  herein  uses  (1)  and  (3)  as  a  basis  of  the  numerical  approximation. 

We  discretize  Q.  with  a  logically  rectangular  mesh  of  hexahedral  cells.  Each  cell  Q  is  the  image  under 
a  trilinear  mapping  of  a  regular  reference  hexahedron  (unit  cube),  Q\  Q  is  determined  by  its  vertices  at 
Vooo,  Vioo,  Voio,  Vooi,  Vno,  Viol,  Von  and  vm,  where  vapy  =  (xapyDaPy,  Zapy)  (Figure  1).  Note  that  the  faces 
of  Q  may  be  non-planar.  The  mapping  associates  with  any  reference  location  ?  =  (x.y.z)  in  Q,  the  point 
r=  (x-.y-.z)  in  Q: 

yax  +  \by  +  xcz  +  \dxy  +  \exz  +  \/yz  +  \gxyz,  (4) 

where  v0  =  v00o,  xa  =  vioo  -  v0,  \b  =  v0i0  -  v0,  vc  =  v00i  -  v0,  \d  =  vno  -  v0  -  va  -  \b,  xe  =  vioi  -  v0  - 
vfl  -  vc,  v/  =  von  -  v0  -  yb  -  vc,  v„  =  vm  -  v„  -  \a  -  xb  -  vr  -  \d  -  -  vf.  This  mapping  extends  the 

conventional  bilinear  mapping  for  a  logically  rectangular  grid  [1,  5],  Note  that  a  fixed  x  in  Q  determines 
a  face  normal  to  the  x  direction,  which  (4)  maps  into  a  corresponding  face  within  Q.  Covariant  vectors, 
defined  as  X(y,£)  =  dr/dx,  Y (x,z)  =  dr/dy  and  Z (x,y)  =  dr/d£,  describe  the  geometry  of  Q.  The  volumetric 
Jacobian  J  for  passing  from  Q  to  Q  is  simply  [6] 

J(x,y,z)  =  X(y,z)  ■  (Y (x,z)  x  Z(x,y)),  (5) 

while,  for  a  face  in  the  logical  x  direction,  the  surface  Jacobian  Tx  becomes  [6] 

F X(x,y,z)  =  |Y(x,£)  x Z(x,y)|,  (6) 

and  the  unit  normal  to  this  surface  is 

n.r  =  (Y  x  Z)  /  |Y  x  Z|.  (7) 

Surface  Jacobians  and  unit  normals  are  similarly  defined  for  faces  in  the  logical  y  and  z  directions.  These 
relations  allow  us  to  define  the  necessary  quantities  used  within  the  CVMFE  method. 
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Figure  2.  Midsection  through  two  adjoining  cells  in  the  logical  x  direction.  Filled  and  open  circles  represent  locations  on  faces 
where  fluxes  are  estimated;  open  circles  represent  flux  locations  on  faces  in  the  logical  y  direction  projected  onto  midsection.  The 
x  symbol  is  used  to  denote  pressure  locations  at  logical  cell  centers. 


3.  CVMFE  METHOD 


The  CVMFE  method  uses  velocity  basis  vectors  Vi.j.tm.n.p  (r)  to  approximate  the  Darcy  velocity  q  in  a  cell 
Qi.j.k  from  the  bulk  cell-face  fluxes.  Here,  i.  j.  k  index  the  cell's  position  in  the  mesh,  while  m.n.  p  determine 
a  face  Fi+m.  j+n.k+P  of  Qi.j.k-  Then  the  approximation  of  q  in  Qi.j.k  can  be  written 


q  ~  fi+l/2J,kVi,j,k-,l/2,0,0  +  fi-l/2J,kVi,j,k,-Q2,0,0  +  fi,j+l/2,kVi,j,k;0,l/2,0 

+  fi,j-l/2,  kVi,j,kfl, -1/2, 0  +  fi,j,k+Q2ViJ,k;0, 0,1/2  +  fi,j,k-\/2 V;j,fc;0,0,-l/2, 


(8) 


where,  for  example,  jG+i/2j/c  is  the  bulk  flux  through  face  Fi+\/2.j.k-  Here,  note  that  “flux”  is  defined  to 
be  the  total  volumetric  discharge  through  a  cell  face,  as  opposed  to  the  volumetric  discharge  per  unit  area, 
which  is  typical  in  the  fluid  mechanics  literature.  Substitute  (8)  into  (1)  and  integrate  both  sides  over  Qi.j.k', 
the  desired  result  is 

fi+\/2,j,k  fi—  l/2,j,k  d"  fi,j+l/2,k  fi ,j—  1  /2,k  "h  fi ,j,k+ 1/2  fi,j,k— 1/2  I  W (x,y,  z)  dxdydz-  (9) 

To  obtain  (9),  by  the  Gauss  divergence  theorem  the  shape  functions  must  satisfy 

/  ^  "Vf  j  km,n,p  dxdyd-Z  —  il.  (10) 


Thus,  the  shape  functions  should  yield  a  discrete  form  (9)  of  the  continuity  equation  (1)  in  Q. 

To  develop  a  discrete  form  for  the  Darcy  relation  (3),  a  control  volume  which  straddles  a  cell  face 
Fi+m.j+n.k+p  is  used;  control  volume  Q*+lp  j  k  for  face  Fi+ij2j^  is  depicted  in  Figure  2.  Weighted  with  a  test 
function  w(+1  /2.j.b  an<l  with  fl  approximated  by  (8),  (3)  is  integrated  over  0*+|.,7  .  k .  The  choice  of  the  test 
function  is  motivated  by  the  efficiency  and  convenience  of  eliminating  face  pressures  adjoining  two  cells; 
see  [4].  Here,  we  use  the  form  proposed  by  Garanzha  and  Konshin  [5]  for  the  test  function,  a  modification 
of  that  originally  described  in  [1],  as  we  have  found  that  it  performs  well  in  most  cases  [4].  Referring  to  face 
Fi+ 1 /2j.k  of  Figure  2,  this  test  function  is 

(  XLJ.k(y.,z)/Ji.J.k(x.,y.,z),  on  Q*+l/4jr, 

w/+ 1  /2,j,k  =  S  ^-i+lJ,k/Ji+l,j,k-i  0,1  Q*i+3/4..j.k- 

( 0,  otherwise, 

where  is  the  covariant  vector  for  dr/dx.  Similar  forms  exist  for  the  logical  y  and  z  directions. 
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Over  the  straddling  control  volume  Q*+l/i  ;  k  the  specific  discharge  approximation  analogous  to  (8) 
involves  the  two  adjoining  cells;  allowing  Yi+l/2j^  to  be  this  approximator,  then 

yi+V2J,k  =  fi+l/2,j,k  1/2, 0,0  +  ^i+lJ,k;-l/2,0,o)  +  fi-l/2J,k^i,j, *;-l/2,0,0 

+  fi+3/2,j,kVi-nj,k;l/2,0,0  +  fi,j+l/2,kVi,j,k;0, 1/2,0  + 

+  fi,j,k+ 1/2  Vi,j,k;  0,0, 1/2  +  fi,j,k- 1/2  Vfj-,*;  0,0, -1/2  +  >/+  1,J+ 1  /2,*  Vi+ 1 1  /2,0 
+  fi+l,j-l/2,kVi+l,j,k;0,- 1/2,0  +/i’+lj,i+l/2V;+ijjit;0, 0,1/2  +  /i+1  ,j,k- 1/2  V,-+ij^;o,0,-l/2-  (12) 

The  discrete  Darcy  relation  is  obtained  from  (3)  by  integrating  (12)  against  (11)  over  £?*+1/o  jk: 

Pi+\,j,k~  Pi, j,k  =  ~P  I  {K-lVi+i/2j,k)-vyi+i/2j,kdxdydz.  (13) 

J  1/2,  j,k 

The  integrations  in  ( 13)  result  in  a  set  of  coefficients  relating  the  bulk  fluxes  on  the  faces  of  Qi.j.k  and  Q/+ 1  ,  j,k 
to  the  difference  in  pressures  at  the  cell  centers  (at  x  —  y  —  z  =  1/2  in  each  cell).  This  method  is  applied  to 
all  interior  faces  of  the  domain,  forming  a  discrete  Darcy  relation  (13)  for  each  face.  The  discrete  continuity 
and  Darcy  equations,  (9)  and  (13),  applied  to  every  cell  within  the  domain,  form  the  discretizations  for  the 
CVMFE  method.  With  irregular  grids,  this  system  of  equations  can  be  non-symmetric,  though  with  certain 
numerical  integration  rules  symmetry  can  be  assured  [5] ;  otherwise,  it  is  similar  in  structure  to  equations 
that  result  from  other  mixed-method  techniques. 

On  a  uniform  grid  with  a  scalar  conductivity  K,  with  the  shape  functions  of  the  usual  mixed  method  and 
with  control-volume  test  functions  as  in  (11),  the  integrations  in  (13)  lead  to  a  tridiagonal  mass  matrix  with 
coefficients  in  the  proportion  of  1/8,  6/8,  1/8.  For  the  usual  mixed  method,  in  which  the  test  functions  are 
the  same  as  the  shape  functions,  the  analogous  proportions  are  1/6,  4/6,  1/6,  and  for  the  usual  finite-volume 
or  finite-difference  methods,  they  arc  0,  1,  0  (diagonal  or  lumped  mass  matrix).  As  pointed  out  in  [5],  based 
on  one-sided  compact  difference  schemes,  the  off-diagonal  coefficient  1/8  of  CVMFE  can  be  expected 
to  preserve  second-order  accuracy  on  non-uniform  grids,  in  which  there  is  no  particular  relation  between 
the  grids  on  the  two  sides  of  an  interface.  This  may  be  the  reason  why  CVMFE  has  calculated  second- 
order-accurate  fluxes  on  non-smooth  quadrilateral  grids  in  2-D,  without  any  need  for  Lagrange  multipliers 
(corresponding  to  face  pressures)  as  in  the  usual  mixed  methods  (e.g.,  [7]). 


4.  FLUX  INTERPOLATION 


This  section  discusses  interpolation  of  cell-face  fluxes  over  cell  interiors,  a  prelude  to  developing  cell- 
velocity  interpolation  functions.  Consider  the  flux  fx(x)  across  an  interior  cell  surface  defined  by  a  fixed 
value  of  x  within  Q:  if  fx  is  the  result  of  a  velocity  field  q(r),  then  this  flux  is: 

fx(x)  =  /q(r)-n xdS  =  [  [  q(r)  •  n{.|Y  x  Z| dydz  =  [  [  q(x,y,z)  ■  (Y  x  Z)(x,y,z)dydz,  (14) 
J X  Jo  Jo  Jo  Jo 

where  nj  is  the  unit  normal  to  this  surface  (7).  The  cross  product  Y  x  Z  is  a  function  of  P  =  ( x,y,z )  and 
represents  a  pole  normal  to  the  surface  at  x .  Its  dependency  on  x  can  be  made  explicit  (see  Appendix  A): 

Y  x  Z  =  ( 1  -  x)  (Yo.y  x  Z0x)  +  x(Y ]x  x  Zix)  -x(  1  -  x)  (Y2x  x  Z2x) ,  (15) 

where  Y0x(£)  =  Y|^=0,  Zttv(y)  =  Z|f=0,  Y Lt(£)  =  Y(f=i,  Zlx(y)  =  Z|*=1,  Y2x  =  Yix  —  Y0.v,  and  Z2x  =  ZU- 

Zox.  Substituting  (14)  into  (15),  the  bulk  flux  across  the  surface  at  x  is 

fx{x)  =  (1  x)fox(x)  +  xf\x(x)  x(l  xjflx^x),  (16) 
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where 


fox(x) 

=  f  f  q (x,y,z)  •  (Y0x  x  Z0x) dydz , 

Jo  Jo 

(17a) 

f\x{x) 

=  [l  [l  q (x,y,z)  •  (Ylx  x  ZLv) dydz, 

Jo  Jo 

(17b) 

f2x(x) 

=  C  C  q (x,y,z)  •  (Y2x  x  Z2.v)  dydz, 

Jo  Jo 

(17c) 

For  a  uniform  flow  field  q(r)  =  q,  ( 17a)-(  17c)  are  independent  of  x.  Since  /o.T  and  f\x  involve  both  Y  and 
Z  at  x  =  0  and  x  =  1 ,  respectively,  and  fox  is  balanced  between x  =  0  and  x  =  1 ,  it  is  natural  to  consider  the 
quadratic  approximation 


QW)  =  (  I  -t)fox(O)  +*/i*(l)  -x(l  -x)/lv(l/2), 


which  is  exact  for  uniform  q.  We  will  also  make  use  of  the  linear  approximation 


LfAx)  =  (l-*)/0x(0)+x/lx(l). 


(18) 


(19) 


The  errors  of  these  approximations  are  derived  in  Appendix  B,  which  we  summarize  here.  To  see  the 
dependence  on  grid  size,  scale  the  x-dimension  to  Ax,  keeping  y  and  z  on  the  unit  interval.  (We  want  to  see 
how  the  flux  error  across  a  fixed  x-face  varies  with  Ax;  refining  it  in  y  and  z  would  have  no  effect,  as  we 
would  integrate  in  y  and  z  over  each  subface  and  would  sum  over  them.)  For  fixed  x,  0  <  x  <  Ax,  the  Green’s 
function  of  the  operator  d2 /dx2  with  homogeneous  Dirichlet  boundary  conditions  is 


G(xX) 


*<Z,<Ax- 


(20) 


Now  f x  —  Lfx  satisfies  these  boundary  conditions  and  ( Lfx )"  =  0,  so  that 

fx(x)  -  Lfx(x)  =  J*G(x&)  (£2  ifx-Lfx)j  (k)dk  =  j\{x£)f"{k)dk-  (21) 

In  Appendix  B,  it  is  shown  that  (21)  leads  to 


{fx~Qfx){x) 

r/xx  ft  I  2  3q 


fj'A  ■  (y'-  x  zi'  -  y»'  x  *-> 


(i  -  •  (Yo*(z)  X  Zftv(y))  +  •  (Yu(z)  x  Zlx(y)) 

h  (qfew)  ~ q  + s  (I  ~  5)  I <Ui)  ~  h  (‘  - 1) 

•  (Y1,(:)  x  (Z2v(v))  1  dydzdt, 


(22) 


and 


r/xx  r  1  r 

(fx  -  Lfx )  (x)  =  (fx  -  Qfx )  (x)  +  /  G(x.X)  / 

J  0  Jo  Jo 

•(Y2a(£)  x  Z2x(y))  dydzdt,. 


1  2  /Ax 


(23) 
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Assuming  smooth  variation  of  q,  Y,  and  Z  with  a,  the  errors  in  (22)  and  (23)  arc  both  0(Ax2),  with 
the  two  powers  of  Ax  coming  from  the  d'q  integral  and  the  Green’s  function.  The  dq/dx  and  q  terms  have 
coefficients  that  are  0(Ax  ')  and  0(Ax  2),  respectively,  but  in  each  case  they  are  multiplied  by  a  difference 
or  product  of  differences  that  restores  the  power(s)  of  Aa. 

Note  that  q  appeal's  in  fx  —  Lfx,  while  fx  —  Qfx  has  only  derivatives  and  differences  of  q.  Thus,  Qfx 
yields  exact  fluxes  for  uniform  q,  while  Lfx  does  not,  though  both  are  second-order.  The  second-order  terms 
in  fx  —  Qfx  are  the  first  two  d2q/dx2  terms,  present  even  on  a  uniform  cartesian  grid  (i.e.,  unavoidable  unless 
q  is  lineal'),  and  the  first  dq/dx  term,  due  only  to  grid  distortion  and  vanishing  if  a- faces  were  uniform.  In 
fx  —  Lfx  there  is  also  the  second-order  term  in  (23). 

The  above  analysis  tacitly  assumes  that  the  coarse-cell  refinement  is  smooth,  e.g.,  along  the  trilinear 
coordinates  set  up  by  the  reference  mapping  (the  coarse  grid  may  be  rough;  only  the  asymptotic  refinement  of 
each  cell  need  be  smooth).  If  so,  then  Y2,  Zi,  and  the  difference  of  Y  x  Z  in  (22)  are  all  O(Ax).  Equivalently, 
in  a  refined  coarse  cell,  the  angular  deviation  from  a  parallelepiped  is  assumed  to  decrease  as  0( Ax)  (so 
that  the  deviation  of  the  vertices  decreases  as  0(Ax2)).  If  not,  i.e.,  if  angular  deviations  do  not  decrease 
(“random”  refinement),  then  Y2,  Zt,  and  the  difference  of  Y  x  Z  are  0(1).  In  fx  —  Qfx ,  the  dq/dx  and  q 
(difference)  terms  become  first-order;  in  fx  —  Lfx ,  the  q  term  becomes  zeroth-order.  Thus,  along  with  its 
exact  interpolation  of  uniform-flow  fluxes  on  distorted  grids,  the  quadratic  flux  interpolation  may  be  crucial 
in  obtaining  convergent  flux  approximations  on  randomly  refined  grids. 

The  next  section  shows  that  the  Piola  transformation,  commonly  used  to  relate  vector-valued  functions 
on  reference  cells  to  irregular  hexahedral  cells,  leads  to  the  linear  flux  interpolation  Lfx.  The  possible 
shortcomings  of  Lfx  noted  here  will  motivate  consideration  of  shape  functions  other  than  those  obtained 
from  the  Piola  transformation. 

In  (18),  the  quadratic  flux  interpolation  requires  that  /o.T(0),  /1x(1)  and  /2*(l/2)  be  known.  Fluxes  /o.r(0) 
and  /i.r(l)  are  components  of  the  solution  vector  for  a  mixed  method.  If  the  cross  product  Ya,-  x  Z2*  in 
fix  is  taken  to  be  a  pole  normal  to  a  surface,  then  |Y2*  x  Zjfdydz  can  be  viewed  as  an  elemental  area  on 
this  surface;  we  adopt  this  and  define  the  surface  as  the  secondary  face,  as  opposed  to  the  primary  faces 
formed  by  the  cell  surfaces.  An  approximation  of  the  secondary  flux  fix{  1/2)  in  terms  of  primary  fluxes  will 
be  presented  later;  here,  it  is  treated  as  known.  If  cells  are  parallelepipeds,  the  last  term  in  (18)  is  null  as 
Y2.V  x  Zix  =  0;  thus,  flux  interpolation  for  regular  grids  will  be  linear,  as  for  the  Piola  transformation. 

In  (17c),  ( y,z )  ranges  over  the  square  (0, 1)  x  (0, 1).  The  secondary  surface  can  be  planar  (vector  Y2_t  x 
Zix  points  in  the  same  direction  for  all  (y,z)),  non-planar  (the  direction  varies  with  (y,£)),  or  twisted  planar 
(the  vector  points  in  opposite  directions  at  ( y,z )  =  (0,0)  and  (1, 1)).  For  twisted  planar  faces,  expressions 
such  as  |  Y2.V  x  Z2x|  dydz  become  suspect  because  such  a  face  would  be  excluded  in  the  2-D  mapping  from  a 
distorted  but  convex  quadrilateral  to  the  unit  square;  such  mappings  are  implicit  to  this  work.  This  behavior 
affects  the  accuracy  of  the  cell-velocity  interpolation  function,  discussed  in  subsequent  sections. 


5.  PIOLA  TRANSFORMATION 

The  Piola  transformation  relates  vector- valued  functions  on  a  hexahedron  Q  to  ones  on  the  reference  cube  Q 
in  a  manner  that  preserves  normal  fluxes  across  primary  cell  faces  and  parts  of  faces.  For  velocity  interpo¬ 
lation  on  hexahedral  meshes,  we  must  define  such  functions  on  Q,  given  certain  values  or  functionals  (e.g., 
integrated  fluxes)  of  these  functions.  Let  r(x,y,£)  =  r(r)  =  (a ,y,z)  be  as  in  (4).  The  Piola  transformation 
takes  a  vector- valued  function  v  on  Q  to  v  on  Q  defined  by 

Drfr) 

v(r(f))  =  7^v(r),  (24) 

where  Dr  is  the  Jacobian  matrix  of  r  and  the  Jacobian  J  is  its  determinant.  Then  [8,  9] 

[  \-n{s)z{s)ds  =  f  x-h(s)z{r{s))ds,  (25) 

JdQ  JdQ 
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where  dQ  is  the  boundary  of  Q.  n(.v)  is  the  outward  unit  normal  vector  at  s  €  dQ,  z.  is  a  scalar- valued  function 
(not  necessarily  continuous)  on  dQ,  and  n(s)  is  the  outward  unit  normal  vector  at  s  €  dQ.  It  is  assumed  that 
v  and  z  are  smooth  enough  so  that  the  integrals  in  (25)  make  sense;  z  can  be  chosen  to  be  1  on  some  parts  of 
dQ  and  0  elsewhere.  Thus,  (24)  preserves  normal  fluxes  across  primary  cell  faces  and  parts  of  faces. 

We  need  to  generalize  (25)  to  interior  faces,  such  as  those  for  a  fixed  x.  The  standard  basis  functions  on 
the  reference  element  Q,  associated  with  primary  face  fluxes,  arc  the  lowest-order  Raviart-Thomas  (RTo) 
functions  [10].  Corresponding  to  the  faces  x  =  0, 1,  y  =  0, 1 ,  £  =  0, 1,  respectively,  arc  the  velocity  fields 
V.v0  =  (1  Vvl  =  (x,0,0),  Vvo  =  (0, 1  —y, 0),  Vvl  =  (0,y,0),  Vz0  =  (0,0, 1  -£),  Vd  =  (0,0,z).  Thus, 

VAi  has  flux  1  across  the  face  Fx\  at x  =  1,  flux  0  across  the  other  five  primary  faces,  and  flux  c  across  interior 
face  Fxc  at  x  =  c.  The  normal  velocity  component  is  constant  on  such  faces. 

Let  v  =  Vri,  and  define  v  by  (24).  Appendix  C  shows  that 

/  vn zds  =  /  cz{r(s))ds,  /  v -nzds  =  0,  (26) 

J  Fxc  •'  J  Fyc 

so  that  fluxes  of  v  across  interior  faces  and  parts  of  such  faces  are  preserved.  By  linear  combinations,  this 
generalizes  at  once  to  any  image  of  an  RTo  function  (Piola-transformed  velocity)  on  Q.  It  follows  that 
the  transformed  normal  component  v  •  n(.v)  on  Fxc  must  be  v  •  n(s)  /Txc{s)  =  c/Txc{s),  pointwise,  where 
Txc  denotes  the  surface  Jacobian  on  Fxc,  so  that  v  •  n  is  proportional  to  1  /Txc.  Consider  the  following  two 
properties  of  a  vector- valued  function  q  on  Q : 

(PI)  The  fluxes  fx{x),  fy(y),  f-(z)  vary  linearly  w  ith  x,  y,  z,  respectively. 

(P2)  The  normal  components  q  n;,  q  n,,  q  •  nz  are  inversely  proportional  to  the  surface  Jacobians  |Y  x  Z|, 
|Z  x  X|,  |X  x  Y|  ,  respectively. 

Then  q  on  Q  is  the  Piola  transform  of  an  RTo  function  q  on  Q  if  and  only  if  q  satisfies  (PI )  and  (P2). 

The  “only  if”  paid  follows  from  the  linear  variation  of  fluxes  of  RTo  functions  on  Q  and  the  flux- 
preservation  properties  (26)  and  generalizations.  For  the  “if”  paid,  let  q  satisfy  (PI)  and  (P2).  On  Q.  there 
is  a  unique  RT0  function  q  with  fluxes  fx( 0),  fx{  1),  /v(0),  fy{  1),  /-( 0),  fz{  1)  on  the  faces  Fx 0,  FxU  Fy0, 
Fy\ ,  Fzo,  Fz i,  respectively.  The  Piola  transform  of  q  satisfies  (PI)  and  (P2),  with  the  same  fluxes  and  normal 
components  as  q  everywhere  on  Q.  Since  a  vector- valued  function  is  uniquely  determined  by  its  three  normal 
components  at  each  point,  q  must  be  the  Piola  transform  of  q.  Note  that  one  consequence  of  this  development 
is  that  the  Piola  transformation  dictates  the  linear  flux  interpolation  in  (19). 

In  2-D,  it  follows  that  for  any  quadrilateral  cell  Q,  a  uniform  flow  q  is  the  exact  Piola  transform  of  an 
RTo  function.  To  see  this,  let  the  vertices  of  Q  be  voo,  vio,  voi,  Vn,  with  the  bilinear  mapping 

r  =  ( 1  —  x)  ( 1  —  y)  voo  +  x(  1  -  y)  vi0  +  ( 1  -  x)yv0i  +xyvn.  (27) 

One  checks  easily  that  X  =  a  +  vb  and  Y  =  c  +  .xd,  for  constant  vectors  a,  b,  c,  d.  The  interior  edge  at  fixed 
x  is  parallel  to  Y  =  (dx/dy,  dy/dy)T,  and  the  flux  of  q  across  this  edge  is  the  integral  of  q  •  n_;]Y|  from  0  to  1, 
where  nA  is  the  unit  normal  to  Y,  and  |Y|  is  the  edge  Jacobian.  Because  n{|Y|  =  (dy/dy.  —dx/dy )T  =  e+xf, 
for  constant  vectors  e,  f,  the  flux  across  the  fixed-x  edge  is  q  •  e  +  xq  •  f,  which  varies  linearly  with  x,  verifying 
(PI).  The  edge  Jacobian  |Y|  is  constant  for  fixed  x,  as  is  the  normal  component  q  •  nA,  verifying  (P2).  Thus, 
in  2-D,  the  linear  flux  interpolation  Lfx  in  (19)  is  exact  for  uniform  flow. 

The  3-D  situation  is  quite  different,  as  shown  by  the  quadratically  varying  fluxes  of  uniform  flow  in  (18). 
A  simple  example  is  the  “truncated  pyramid”  hexahedron  Q  (Figure  3a),  whose  bottom  face  atz=  {sq  —  s\)/2 
(z  =  0)  is  a  square  of  side  so  (.—sq/2  <x<  sq/2,  —sq/2  <  y  <  sq/2),  and  whose  top  face  at  z  =  (si  —  sf)/2 
(z  =  1)  is  a  squai'e  of  side  s\  (—s\/2  <x<  s\/2,  —s\/2  <  y  <  sq/2)\ 

x  =  a(2x  —  l)/2,  y  =  a( 2y-l)/2,z  =  (z- l/2)(si -s0),  (28) 

where  a  =  .so  +  —  so).  The  four  lateral  faces  are  planar  trapezoids  (see  Figure  3a).  For  vertical  uniform 

flow  q  =  (0,0, 1),  the  flux  across  horizontal  interior  faces  for  fixed  z  is  the  cross-sectional  area, 

fz(z)  =  [(l-z)*>  +  fci]2-  (29) 
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a.  Truncated  pyramid.  b.  Tent  with  curved  roof. 

Figure  3.  a.  Hexahedral  cell  in  form  of  a  truncated  pyramid.  Uniform  flow  parallel  to  vertical  z  axis,  as  indicated  by  arrow,  is 
assumed,  sq  and  indicate  the  length  of  the  lower  and  upper  edges,  respectively,  b.  Hexahedral  cell  in  form  of  a  tent  with  curved 
roof.  Uniform  flow  parallel  to  horizontal  x  axis,  as  indicated  by  arrow,  is  assumed.  Dash-dot  lines  show  maxima  and  minima  of 
saddle  surface  forming  cell  face  z  =  1 .  For  both  figures,  vertices  vooo  and  vi  1 1  are  indicated;  see  Figure  1  for  orientation. 

A  Piola- transformed  RTo  function  with  the  correct  top  and  bottom  fluxes  has  linearly  varying  flux 

Lfz(z)  =  (l~z)M0)  +  zfz(l)  =  (1  -z)sl  +  zs\.  (30) 

As  (29)  and  (30)  do  not  agree,  Piola-transformed  RTo  functions  cannot  represent  a  uniform  flow  field  exactly 
in  3-D,  even  on  such  a  simple  element.  Quadratic  flux  interpolation  (18)  treats  this  case. 

Another  example  is  a  “tent  with  curved  roof’  (Figure  3b)  with  unit  square  base  at  z  =  0  (£  =  0),  vertical 
lateral  faces,  height  ho  at  two  opposite  vertices  (0,0,/zo),  (1,  I .  ho ) ,  and  height  h\  at  the  other  two  vertices: 

x  =  x,  y  =  y,  z  =  {ho[(l-x){l-y)+xy]+hi[{l-x)y  +  x(l-y)])z.  (31) 

For  horizontal  uniform  flow  q  =  (1,0,0),  by  symmetry  the  flux  across  the  curved  roof  Fz\  at  £  =  1  is  0 
(this  also  follows  from  the  divergence  theorem).  An  arbitrary  Piola-transformed  RTo  function  v  with  flux 
0  across  F-\  is  a  linear  combination  of  the  transforms  VYo,  Vxj ,  Vyo,  Vv i ,  V-o,  Vzi  of  the  six  reference  basis 
functions.  By  (25),  the  transforms  other  than  Vzi  have  vanishing  normal  component  on  Fz\,  pointwise.  Thus, 
the  contribution  of  V-i  to  v  has  flux  0  across  Fz\,  so  that  its  coefficient  in  the  linear  combination  is  0;  hence 
v  •  n  vanishes  on  Fz i,  pointwise.  Since  q  •  n  ^  0,  q  is  not  a  Piola-transformed  RTo  function.  More  generally, 
a  3-D  shape-function  space  with  (a)  degrees  of  freedom  corresponding  to  face  fluxes  and  (b)  pointwise 
preservation  of  null  normal  components  across  surfaces  cannot,  for  general  hexahedra  that  arc  tril i near 
images  of  a  reference  cube,  contain  the  constant  vector  fields.  This  example  points  out  issues  of  pointwise 
velocity  interpolation,  as  opposed  to  integrated  flux  interpolation,  and  the  difficulties  of  treating  hexahedra 
with  non-planar  faces. 


6.  VELOCITY  INTERPOLATION 

Given  primary  face  fluxes  /o.T(0),  etc.,  and  secondary  fluxes  /2*(l/2),  etc.,  this  section  determines  interpo¬ 
lated  velocity  approximations  at  any  point  r  in  a  cell  Q.  We  build  these  approximations  from  components 
that  parallel  the  covariant  vectors  X,  Y  and  Z.  For  the  exact  Darcy  velocity  q(r)  these  components  arc 
designated  qx(f),  qy(r)  and  qz(r),  with 

q  =  Qx  +  qy  +  qz-  (32) 
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The  magnitude  Ux( r)  of  the  velocity  normal  to  an  arbitrary  surface  at x,  with  unit  normal  nv,  is 


US  =  q  •  n f  =  qx  nf, 


(33) 


because  nA  •  qy  =  n?  •  qz  =  0-  Conversely,  if  the  normal  component  Uz(x.y.z)  is  specified,  then 


„  Us(x,y,z)\YxZ\(x,y,z)„,^,  „X|Vwry| 

qx  {x,y,z)  =  ,Vv7U<.,'  X(y,z)  =  Ux-  Y  x  Z 

X(y,z)  •  (Y  x  Z)(x,y,z)  J 


(34) 


satisfies  (33).  Thus,  in  theory,  a  shape  function  can  match  an  arbitrary  U$,  and  similarly  for  qy  and  qz- 
Equation  (34)  is  too  general;  practical  shape  functions  that  specify  Ux{x,y,z)  must  be  defined.  The  rel¬ 
evant  unknowns  for  CVMFE  are  primary  face  fluxes  /o.T(0),  fix(  I ),  /ov(0),  /iy(l),  /o-(0),  f\ - ( I ) ;  hence, 
Ux(x,y,z )  should  depend  on  fox{ 0),  f\x{  I ),  and  possibly  the  others.  The  Piola  transformation  takes 


TJ  r  .  (l-x)fox(0)+xMl)  Lfx(x) 

*  'yv  J  |Y  x  Z\(x,y,z)  |Y  x  Z| {x,y,z)  ’ 

which  we  verify  by  checking  properties  (PI)  and  (P2).  For  (PI),  as  in  (14), 


(35) 


fx(x)=  ['  [  U£{x,y,z)\YxZ\{x,y,z)dydz  =  Lfx(x),  (36) 

Jo  Jo 

and  (P2)  is  immediate.  The  potential  shortcomings  of  the  Piola-transformed  shape  functions  have  been  noted 
in  previous  sections,  and  we  therefore  consider  alternatives  to  properties  (PI)  and  (P2). 

The  fluxes  of  uniform  flow  can  be  matched  exactly  if  we  replace  (PI)  with 


(Ul)  The  fluxes  fx(x),  fy(y),  fz(z )  vary  quadratically  with  x,  y,  z,  respectively,  in  the  manner  of  (18). 

Here  f2x(  I  /2)  must  be  determined  by  /oA(0),  f\x(  I ),  and  possibly  the  other  primary  fluxes,  in  ways  to  be 
described  later.  Because  a  uniform  flow  has  constant  normal  component  on  a  planar  face,  consider  also 


(U2)  The  normal  components  q  •  n^,  q  •  nA,  q  •  m  arc  constant  on  interior  faces  with  fixed  x,  y,  z,  respectively. 
For  a  fixed-x  interior  face,  with  flux  fx(x),  (U2)  requires  that 


Ui{x,y,z)=fx{x)/Ax{x)  (37) 

(unlike  (P2),  for  which  Ux(x,y,z)  =  fx(x) / 1 Y  x  Z|  (Jc,y,z)),  where  Ax  is  the  surface  area  of  the  interior  face: 

Ax{x)  =  [l  [l  |Y  x  Z| (x,y,z) dydt  (38) 

Jo  Jo 

Only  if  the  interior  face  is  planar  is  this  straightforward  to  compute  for  all  x,  making  (U2)  somewhat  im¬ 
practical.  It  also  matches  uniform  flow  pointwise  only  under  those  circumstances.  (In  general,  if  the  normal 
direction  changes  by  0(1)  across  a  hexahedral  cell  of  size  0(1),  as  it  can,  then  upon  tril i near  refinement  the 
direction  will  change  by  0(Ax)  across  a  refined  hexahedron  of  size  O(Ax);  thus,  first-order  accuracy  in  the 
pointwise  representation  of  uniform  flow  is  the  best  that  one  can  hope  for.)  Accordingly,  we  formulate  the 
less-restrictive  property  that  requires  constancy  at  the  ends  (primary  faces)  only: 


(UE2)  The  normal  components  q  •  nj,  q  •  n-p,  q  •  m  arc  constant  on  the  primary  cell  faces. 


Property  (UE2)  does  not  specify  how  normal  components  vary  in  the  interior;  this  variation  will  be  a  by¬ 
product  of  combinations  of  (UE2)  with  the  flux  variations  (PI)  and  (Ul). 

A  simple  modification  of  the  Piola  shape  functions  that  suffices  to  represent  uniform  flow  on  the  “trun¬ 
cated  pyramid”  is  the  replacement  of  linear  with  quadratic  flux  interpolation,  combining  (Ul)  with  (P2): 


Ux(x,y,z)  = 


(1  —  .y)/o.t(0)  ±xMl)  -Jf(l  -x)f2x(l/2) 

|YxZ|  (x,y,z) 


Qfx(x) 


|Y  x  Z\(x,y,z) 


(39) 
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For  elements  with  variable  surface  Jacobians  on  primary  and/or  interior  faces,  combine  (PI)  with  (UE2): 


Ux{x,y,z)  =  (1-4) 


/ftv(O) 

4,(0) 


0x+Ax(\)  ]x: 


(40) 


where 


OCftr 


Yq.v  X  Z0x\(y,z) 
|YXZ|  (x,y,z)  ’ 


«1,= 


Yu  x  Zix\(y,z) 
|YxZ|  (x,y,z)  ' 


(41) 


At  4  =  0  and  4=1,  (40)  reduces  to  the  form  of  (37),  with  constant  normal  component.  The  interpolated 
flux  is  as  in  (36),  which  is  the  reason  for  the  factors  (Xqx  and  ai,  in  (40).  An  analogous  development  with 
quadratically  varying  flux  combines  (Ul)  with  (UE2)  to  obtain 


j,  ~\/o,(0)  . 

Ux{X,y,Z)  —  (1  x)  —  -7—r-  OtOx  +  X—  7—rOCl_v  X{1  Xj  -  OC2.V, 


4,(0) 


4,(1) 


42, 


(42) 


where  a 2,  and  A2x  are  discussed  below.  For  the  secondary  surface,  two  alternatives  are  suggested  for  0C2,: 


_  |Y2,xZ1v|  _  (Y2,xZ2,)-W2, 

a2x |Y  x  Z|  ’  azr_2  |W2,||YxZ| 


(43) 


where  Y2,(£)  =  Yl,  -  Y0„  Z2,(y)  =  Zlx  -  Z0„  and  W2,  is  defined  as 

W2,=  [ ‘  [\y1x(z)  xZ2x(y))  dydz  =  C  Y  2x{z)dzx  f'z^dy  =  Y2,(l/2)  x  Z2,(l/2).  (44) 

Jo  Jo  Jo  Jo 

To  show  the  second  equality  in  (44),  write  out  each  component  of  the  cross  product  in  the  vector  double 
integral,  observe  that  the  terms  in  the  resulting  scalar  double  integrals  can  be  factored  into  scalar  single 
integrals,  and  finally  reassemble  these  into  the  components  of  a  cross  product  of  vector  single  integrals.  The 
last  equality  follows  from  the  lineality  of  Y2,  and  Z2,  with  respect  to  z  and  y,  respectively.  Related  to  these 
two  alternatives  are  two  alternative  definitions  for  the  area  42,  of  the  secondary  surface: 

42,-1  =  [ 1  [l\Y2xxZ2x\dydz,  42,-2  =  ['  [ ' ( Ylx  x  Z2v)  •  dydz  =  |W2,|.  (45) 

Jo  Jo  Jo  Jo  |W2,| 


Both  expressions  in  (43)  are  approximate  for  non-planar  secondary  faces.  For  twisted  planar  secondary 
faces,  42,_i  is  incorrect  because  it  is  based  on  a  mapping  for  a  convex  region  and  prevents  cancellation 
of  opposite  directions  in  Y2,  x  Z2,;  42,_2  is  correct  in  this  case,  and  a2x  2  should  be  superior  to  0C2,-i  as 
the  twisted-face  differential  area  is  used  in  the  construction  of  this  ratio.  For  planar,  non-twisted,  secondary 
faces,  (Y2,  x  Z2x)(y,z)  is  a  positive  multiple  of  W2,,  so  that  both  forms  in  (43)  and  (45)  produce  identical 
results. 

For  the  velocity  magnitude,  though  other  forms  have  been  studied,  we  have  found  that  (42)  performs  well 
because  the  ratios  (41)  and  (43)  adjust  the  magnitude  by  comparing  a  differential  element  on  the  interior 
surface  with  one  on  the  cell  surface.  In  contrast  to  Ux(x,y,z)  of  (33),  the  approximate  velocity  magnitudes 
in  (40)  and  (42)  become  constant  on  the  cell  faces  at  x  =  0  and  x  =  1.  This  results  because  (40)  and 
(42)  are  based  on  the  flux  interpolation  functions  (19)  and  (18),  respectively;  these  interpolation  functions 
assume  bulk  fluxes  at  the  cell  faces  to  approximate  the  bulk  flux  on  a  congruent  interior  face.  This  level  of 
approximation  is  consistent  with  the  discrete  continuity  and  Darcy  equations,  (9)  and  (13). 

Use  U_x  from  (42)  to  approximate  (34);  then  the  approximation  Vx  of  qx  in  (34)  is 


Vx  =  j  [( 1  -  x)  +  x (3,i/,i  -  x{  1  -  x)  (3 ,2/,2] ,  (46) 

where  /,o  =  /o,(0)  and  fx\  =  /i,(l)  are  the  actual  fluxes  across  cell  faces  x  =  0  and  4=1,  and  fx2  is  the 
secondary  flux  to  be  determined  later.  The  ratios  (3,o  and  (3,i  and  the  suggested  alternatives  for  (3,2  are 


P.v0 


YxZ 


(3,i 


i=0 


YxZ 


x=l 


(3,2  1 


Y2r  x  Z2,| 
42,-1 


[3,2—2 


(Y2,  x  Z2,)  •  W2, 
|W2,|2 


(47) 
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Similar  expressions  can  be  derived  for  interpolation  functions  in  the  Y  and  Z  directions: 

Vy  =  J  [( 1  -  y)  Pyo/yO  +  y  Pvl/yl  -  y ( 1  -  y)  Py2/y2] , 

Vz  =  y[(l-z)Pzo/zO  +  Zpzl/zl-z(l-z)Pz2/z2], 

where 


PyO  - 


Z  X  X 


Ay 


Pyl  - 


y=o 


ZxX 


P;0  - 


y=l 


XX  Y 


Psl  - 


z=0 


Xx  Y 


l=t 


(48) 

(49) 


Py2— 1 


Pz2-1 


Z2y  X  X2y| 

A2y— 1 

|X2;  X  Y2-l 
A2z-  i 


Py2— 2  = 

Pz2-2  = 


( Z2y  X  X2>,)  •  W2>. 
|W2y|2  ’ 

(X2;  X  Y2t)  •  W2; 

|W2z|2- 


The  associated  secondary  vectors  arc  X2>,  =  X|y=i  —  X|y=o,  Z2v  =  Z|y=i  —  Z|-p=o,  X2,  =  X| s=i  —  Xp=o  and 
Y2,  =  Y|*=i  —  Y|*=o;  fyo,  fy i,  fzo  and  fz\  arc  fluxes  across  cell  faces  at  y  =  0.  y  =  1,  z  =  0  and  2=  1;  and 
fy2  and  fz2  are  fluxes  associated  with  the  secondary  faces  in  the  logical  y  and  z  directions.  Areas  Ay(y)  and 
Az(z)  are  similar  to  (38)  and  secondary  areas  4 2v  |  and  A2z_j  arc  similar  to  42x  i  in  (45). 

The  cell  approximation  Vc  of  the  velocity  q  can  be  obtained  from  these  relations  as 


Vc  =  Vx  +  Vy  +  Yz. 


(50) 


For  uniform  flow,  where  (18)  is  exact,  Vc  can  equal  q  provided  that  all  cell  faces  are  planar;  however, 
because  of  possible  errors  in  the  velocity  magnitude  correction  ratios  Px2,  Pv2  and  p,2  for  the  secondary 
face,  this  equality  may  not  hold.  This  discrepancy  will  be  discussed  further  subsequently. 


7.  VELOCITY  SHAPE  FUNCTIONS 

This  section  finalizes  the  shape  functions  by  developing  accurate  approximations  for  the  secondary  fluxes 
fx 2,  fy2  and  f-2  in  terms  of  the  primary  fluxes  fx q,  fx\ ,  fyo,  fy\ ,  f-o  and  f-  \ .  The  exact  relation  for  /x2  is  (17c) 
with  x  =  1/2;  one  could  propose  (50)  to  approximate  q  in  (17c),  except  that  Vc  itself  requires  knowledge 
of  fx2,  fy2  and  f-2.  We  propose,  instead,  to  approximate  q  in  (17c)  with  a  bulk  value  (q)  for  the  cell  Q.  The 
approximation  /x2  for  the  secondary  flux  /v2  would  then  be 

fx 2  =  (q>  •  [ 1  /  WxZ 2x)dydz  =  (q)  •  W2x,  (51) 

Jo  Jo 

where  W2jt,  derived  in  (44),  is  a  bulk  equivalent  of  Y2x  x  Z2v  that  also  appeal's  in  (43)  and  (45)  to  correct 
for  twisted  planar  secondary  faces.  The  fact  that  |W2v|  gives  the  correct  area  for  a  twisted  planar  secondary 
face  inspired  the  differential  area  used  in  these  expressions.  For  the  bulk  cell  velocity  (q),  we  propose: 

(q)  =  ^[Vc(0,0,0)  +  VC(1,0, 0)  +  Vc(0, 1,0)  +  Vc(0,0, 1) 

+  Vc(l,  1, 1)  +  Vc(0, 1, 1)  +  VC(1,0, 1)  +  Vc(l,  1,0)].  (52) 

That  is,  the  cell  velocity  Vc  at  the  corners  is  averaged.  At  the  corners,  the  secondary  fluxes  are  not  needed 
in  Vc,  so  that  we  can  write 

Vc(i,j,k)  =  [fM-.k)  (3  A)  |  fyjY(i.k)  (3W(/,A)  +/z*Z(i,  j)  MiJ)]/J(iJ,k),  (53) 
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where  i  —  0, 1;  j  —  0, 1;  k  —  0, 1.  With  (51),  (52)  should  produce  a  reasonable  estimate  of  the  secondary 
flux,  fx2 ;  for  uniform  flow,  it  is  exact.  Appendix  D  shows  that,  for  a  general  smooth  velocity  field  q,  this 
approximation  preserves  the  second-order  accuracy  of  the  quadratic  flux  interpolation  Qfx,  as  in  (22). 

This  shows  that  fx 2  can  be  approximated  as  a  linear  combination  of  the  primary  fluxes  fx 0,  fx  1,  fy 0,  fy\ , 
fz 0  and  f- 1 ,  and  similarly  for  fy 2  and  fz 2.  The  general  secondary-flux  approximation  is,  for  K  =  x,  y,  or  z,  and 
/  =  0, 1 : 


/k2 


r Kx<)fx0  +  To:  1  fx  I  +  r KvO/vO  +  Tcvl/vl  +  rKz0fz0  +  rK-l/;l  5 


(54) 


Oct/  — 


w2k 


Ocv/  = 


w2k 


_  y  y  M»b/j)X(m,/j) 

_  y  y  ^yl(m,p)Y{m,p) 
=  WlK  _  y  y  $-l{m,p)Z{m,p) 

8  m=0P=0  ■'Kp.O 


(55a) 

(55b) 

(55c) 


Here,  J(x,y,z)  is  the  Jacobian  in  (5)  and  |3V/,  (3V/  and  [3-/  are  as  in  (46)-(49).  Equation  (54)  follows  by 
substituting  (53)  into  (52)  and  (52)  into  (51).  Following  (44),  the  bulk  vectors  W2v  and  W2-  arc 


W2v= 


W2z= 


fQ  J\z 2y(x)  x  X2y(z))dxdz  =  Z2y  (1/2)  x  X2v  (1/2) , 

f 1  I' (X2z(y)  x  Y*(x))  dxdy  =  X* (1/2)  x  ¥*(1/2) . 
Jo  Jo 


(56a) 

(56b) 


These  forms  complete  the  velocity  interpolation  functions  Vx,  Vy  and  V/  in  (46)-(49). 

Now  shape  functions  can  be  obtained  directly  from  (50)  by  collecting  coefficients  of  fx o,  fx\ ,  /vo,  /yi ,  /-o 
and  f- [ .  In  terms  of  shape  functions  v.ro,  VA  i ,  Vvq,  Vyi,  V-o  and  Vzi,  the  cell  velocity  becomes 


Yc  —  fx OYtO  T  fx\  Vtl  T  fytiVyO  T  fy lVyi  +  /jOVjO  T  /zlWlj 


(57) 


where 

V.vo  =  j  [( 1  —  x)X  (3^0  —  x(  1  —  x)  X  $xirxxO  —  y(  1  —  y)  Y  $y2ryx0  —  £(  1  —  £)  Z  (58a) 

V.vi  =  J  [J?x  P,1  -  jf  ( 1  -  x)  X  (3x2rxxl  -  y ( 1  -  y)  Y  $y2ryxl  -  £(  1  -  z)  Z  (3z2rwl ] ,  (5 8b) 

VvO  =  j[(l  — y)Y(3y0  — x(l  — f)  X(3  x2rxyo  y(l  y )  Y  (3v2 ryv0  2(1  z)Zfiz2Fzyo\i  (58c) 

Vv,  =  j[yY(3,t  -x(l -x)Xfarxyi -y{l -y)  Y^i  -£(1  -£)  Zp^],  (58d) 

VjO  =  j[(l— £)Z(3z0-i(l-x)X(3x2rxz0— y(l— y)Y(3y2ryzo-£(l-£)Z(3,2rzz0],  (58e) 

Vzi  =  j[£Zpzl  -x(l  -i)X(3x2rxd  —  y(l  -y)  Y(3y2ryd  -£(1  -£)ZPz2rzd],  (58f) 


Properly  reconfigured,  Vc  from  adjoining  cells  approximate  V,+i/2jd  of  (12).  For  a  cell  Qi.j.k,  the  shape 
functions  arc  as  in  (8)  and  comply  with  property  ( 10) ;  that  is,  a  discrete  continuity  equation  is  formulated  for 
Qi.j.k-  This  completes  the  CVMFE  method.  Note  that  if  normal  components  were  approximated  with  linear 
forms  such  as  (40)  instead  of  (42),  then  only  the  first  term  in  each  of  (58n)-(58/)  would  remain  (also  see 


final_report.tex;  6/11/2001;  9:04;  p.12 


13 


reference  [4]).  Because  of  the  complexity  of  the  shape  functions  and  the  test  function  (11),  the  integrations 
in  (13)  were  calculated  with  a  very  precise  quadrature;  the  effect  of  quadrature  will  be  investigated  in 
the  future.  The  resulting  algebraic  equations  were  solved  with  the  Schur  complement  to  decouple  the  flux 
equations  from  the  pressure  equations,  followed  by  a  variant  of  the  iterative  scheme  of  [11]. 

The  shape  functions  for  3-D  hexahedral  cells  can  have  a  quadratic  dependence  when  the  secondary  fluxes 
arc  nonzero.  This  condition  occurs  whenever  a  cell  has  non-parallel  opposite  faces;  only  regular  meshes  or 
irregular  meshes  with  a  2-D  aspect  will  retain  strictly  linear  shape  functions.  For  uniform  flow,  these  shape 
functions  can  give  exact  velocity  results  if  the  hexahedral  cell  has  planar  primary  (exterior)  and  secondary 
faces.  Mesh  construction  can  usually  produce  planar  primary  faces,  but  planar  secondary  faces  can  be  more 
problematic,  and  the  distortion  on  the  secondary  faces  can  be  rather  severe.  In  particular,  secondary  faces 
can  be  non-planar  even  when  all  exterior  faces  are  planar.  Planar  secondary  faces  can  be  forced  by  insisting 
that  each  cell  have  at  least  one  pair  of  parallel  opposite  faces;  however,  these  secondary  faces  may  be  twisted 
as  well  as  planar. 


8.  UNIFORM  FLOW  TESTS 

Uniform  flow  should  be  approximated  accurately  in  any  numerical  simulation  of  flow  through  porous  media, 
and  CVMFE  usually  does  so.  However,  accuracy  can  be  lost  in  certain  cases  of  irregular  cells,  particularly  if 
either  the  primary  or  secondary  cell  faces  arc  non-planar.  Consider  non-planar  primary,  exterior  faces.  Two 
of  the  three  covariant  vectors  arc  tangent  to  every  exterior  face,  and  the  velocity  shape  functions  depend  on 
these  vectors,  so  that  approximate  velocities  near  non-planar  faces  should  be  curved  along  with  these  faces, 
deviating  from  a  uniform  flow  field.  The  “tent  with  curved  roof”  in  Section  5  is  an  example.  For  secondary 
faces,  the  reasoning  is  less  clear  but  likely  similar. 

The  error  between  simulations  and  exact,  uniform  flow  is  measured  with  the  L?  norm 


where,  for  all  faces  Fj  e  Q.,  J)  and  /,•  are  the  estimated  and  true  fluxes.  A,  is  the  area,  and  Q*  is  the  control 
volume  straddling  Fj.  The  numerator  of  (59)  is  a  discrete  volume  integral  of  the  square  of  the  velocity  error, 
and  the  denominator  normalizes  the  error  to  be  independent  of  the  domain  volume.  In  the  test  problems  with 
a  dominant  linear  flow  direction,  the  mean  velocity  will  be  1;  in  the  next  section,  for  a  non-uniform  flow 
test  with  corner-to-corner  flow,  it  will  be  of  order  1 .  Thus,  the  reported  L2  norms  can  be  regarded  as  relative 
errors  and  can  be  compared  from  one  test  to  another. 

Three  different  coarse  meshes  were  used  in  the  uniform  flow  tests:  in  tests  1  and  2,  cells  were  quasi¬ 
random  with  planar  primary  faces,  while  the  vertices  of  test  3  were  random  and  all  interior  primary  cell  faces 
were  non-planar.  The  domain  fl  was  a  regular  cube  in  all  cases  with  4x4x4  coarse  mesh.  Each  discretization 
was  refined  smoothly  through  three  levels.  At  each  level,  each  cell  was  divided  into  eight  subcells,  halving  it 
in  the  reference  cube  in  each  direction.  This  yielded  refined  meshes  of  8x8x8,  16x  16x  16  and  32x32x32. 

In  test  1 ,  irregular  cells  with  planar,  quasi -random  faces  were  generated,  starting  with  random  vertices 
on  three  adjoining  exterior  faces  of  £2.  These  vertices  were  created  from  a  regular  mesh  on  each  face  by 
multiplying  each  regular  vertex  location  by  a  random  number  close  to  one.  Starting  at  the  corner  where  the 
three  exterior  faces  meet,  cells  were  created  by  connecting  vertices  of  three  existing  cell  faces  to  an  interior 
vertex,  requiring  that  the  new  cell  faces  forming  this  connection  be  planar.  This  results  in  interior  cells  with 
planar,  quasi-random  primary  faces.  The  secondary  faces  of  cells  in  this  grid  are  rarely,  if  ever,  planar.  In  test 
2,  the  random  vertices  on  the  three  adjoining  exterior  faces  of  £2  were  required  to  fall  on  parallel  grid  lines 
perpendicular  to  the  x  direction  (Figure  4).  This  affects  primarily  the  x  cell  faces,  which  arc  all  parallel  and 
follow  the  grid  lines  perpendicular  to  the  x  direction.  The  y  and  z  cell  faces  remain  quasi-random  after  the 
inward  projection  procedure,  but  arc  constrained  in  sharing  vertices  with  the  x  faces.  The  chief  difference 
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Figure  4.  Three  exterior  faces  of  the  test  2  domain  that  share  vertex  (0,0,0).  Fold  up  the  left  (yz )  face  along  the  y-axis  and  the 
bottom  (xz)  face  along  the  v-axis. 


between  the  two  types  of  grids  is  that  the  secondary  faces  of  cells  in  the  test  2  grid  are  always  planar  but 
occasionally  twisted.  The  coarse  mesh  for  test  3  was  created  by  randomly  perturbing  the  vertex  locations 
of  a  regular,  4x4x4  grid;  the  vertices  on  dQ  were  perturbed  only  within  the  plane  forming  the  surface. 
None  of  the  interior  cell  faces  were  planar.  The  random  perturbation  was  limited  to  a  maximum  of  5%  of  the 
unperturbed  cell  dimension;  this  small  amount  of  grid  distortion  may  easily  be  expected  in  flow  simulations. 

In  test  1 ,  where  all  x,  y  and  z  are  quasi-random,  the  uniform  flow  q  =  ( 1 , 0, 0)  was  imposed  via  homo¬ 
geneous  K  and  appropriate  flux  boundary  conditions.  It  is  then  simple  to  calculate  the  exact  fluxes  across 
the  interior  faces  of  the  discretization.  In  test  2,  with  quasi-random  y  and  z  faces,  similar  flux  conditions 
were  imposed  to  obtain  q  =  (1,0,0)  and  (in  a  separate  calculation)  q  =  (0,1,0).  These  two  different  flux 
conditions  were  used  because,  on  the  average,  the  test  2  mesh  is  anisotropic.  The  uniform  flow  in  test  3  was 
diagonal  across  the  mesh,  from  cell  Q\\  \  to  cell  Q  444.  Results  from  the  non-linear  velocity  shape  functions 
(sfl,  sf2)  developed  herein  are  compared  with  linear  velocity  shape  functions  (sfO);  explicit  forms  for  sfO 
can  be  found  in  [4] .  Here  sfl  and  sf2  refer  to  secondary-flux  correction  factors  (3t2  1 ,  (3V2  1 ,  (3-2  1  and  (3X2  2, 
(3v2  2,  Pz2-2-  respectively,  in  (46).  The  L2  norm  results  are  in  Table  I. 

For  test  1 ,  the  non-linear  shape  functions  (sfl ,  sf2)  arc  superior  to  the  linear  ones  (sfO)  over  the  entire 
range  of  mesh  refinement  (1  x-8x)  by  one  to  two  orders  of  magnitude.  Thus,  the  quadratic  flux  interpolation 
(18)  has  a  substantial  impact  on  the  second-order  error  in  velocity  interpolation,  particularly  on  coarse  grids, 
which  arc  important  in  practice.  For  the  non-linear  shape  functions,  there  is  little  difference  between  the  two 
correction  factors,  and  results  become  identical  as  the  mesh  is  refined.  For  either  sfl  or  sf2,  the  L2  norms 
are  not  large  in  any  case;  this  probably  reflects  the  moderate  degree  of  distortion  contained  in  the  mesh. 
With  the  first  and  subsequent  refinements,  the  cell  faces  become  slightly  non-planar,  adding  some  error  to 
the  solution.  This  may  explain  the  slightly  smaller  reduction  factors  for  the  first  refinement.  Subsequent 
refinements  reduce  the  L2  norm  by  a  factor  of  3.3  to  3.7,  as  the  effect  of  non-planarity  decreases.  The 
reduction  factors  for  sfO  are  5  to  5.5;  as  the  secondary  flux  becomes  less  important  with  refinement,  the  cells 
become  closer  to  parallelepipeds,  and  sfO  l 2  norms  come  closer  to  those  of  sfl  and  sf2. 

For  test  2,  with  parallel  faces  in  the  x  direction  and  flow  in  the  x  direction  (x  fd),  the  L 2  norm  results 
for  the  coarse  grid  (lx)  arc  two  orders  of  magnitude  greater  for  sfO  than  for  sfl.  However,  after  three  mesh 
refinements  (8x),  the  ratio  between  sfO  and  sfl  is  only  a  factor  of  three;  the  more  rapid  reduction  of  this 
ratio  here,  relative  to  the  test  1  case,  probably  reflects  the  strongly  diminished  importance  of  the  secondary 
flux  with  mesh  refinement.  For  x  fd,  the  non-linear  shape  functions  sfl  and  sf2  give  very  different  results  on 
the  coarse  grid  (1  x);  sf2  reflects  an  exact  solution.  This  does  not  occur  for  sfl  because  the  random  grid,  in 
this  case,  produces  planar  x-direction  secondary  faces,  some  of  which  are  twisted.  While  the  sf2  correction 
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Table  I.  L r  norm  results  for  uniform  flow  tests.  Test  1:  all  cell  faces  random  and  planar;  flow  in  a  direction. 
Test  2:  cell  faces  y  and  z  random;  x  faces  parallel;  all  cell  faces  planar.  Test  3:  all  vertices  random;  all 
internal  faces  non-planar;  diagonal  flow.  sfO:  linear  shape  function  [4],  sfl:  non-linear  shape  function  using 
first  secondary  flux  correction  factors  $X2-U  etc-  sf2:  same  with  (3f2-2>  etc.  fd:  flow  direction.  *:  <  10-1,) 


test 

grid  refinement,  L2  norm 

no 

case 

lx 

2x 

4x 

8x 

sfO 

3.591  x  nr04 

7.103  x  10“05 

1.286  x  10“05 

2.344  x  10“06 

1 

sfl 

6.098  x  10-06 

2. 1 1 1  x  10“06 

6.453  x  10-°7 

1.737  x  10“07 

sf2 

6.116x  10“06 

2. 1 1 1  x  10“06 

6.452  x  10“07 

1.737  x  10“07 

sfO,  x  fd 

1.817  x  10-03 

3.176  x  10-04 

5.777  x  10“05 

1.070  x  10“05 

sfl ,  x  fd 

1.862  x  10“05 

3.608  x  10“05 

1.202  x  10“05 

3.311  x  10“06 

2 

sf2,  x  fd 

exact* 

3.599  x  10“05 

1.202  x  10“05 

3.311  x  10“06 

sfO,  sfl. 
sf2,  y  fd 

exact* 

1.666  x  10-05 

5.210  x  10-06 

1.401  x  10-06 

sfO 

1.232  x  10“04 

4.045  x  10“05 

1.151  x  10“05 

3.199  x  10-06 

3 

sfl 

1.039  x  10“04 

3.940  x  10“05 

1.137  x  10“05 

3.182x  10“06 

sf2 

1.038  x  10-04 

3.940  x  10“05 

1.137  x  10“05 

3. 182  x  10-06 

factors  can  compensate  completely  for  twisted  planar  faces,  the  sfl  correction  factors  cannot.  Both  non¬ 
linear  shape  functions  produce  nearly  identical  L?  norms  with  the  first  refinement  (2  x )  of  the  mesh,  and 
both  degrade  relative  to  the  initial  mesh  (lx).  This  is  likely  due  to  the  introduction  of  non-planar  primary 
faces  in  the  2x  refined  mesh.  With  refinement,  the  tendency  for  a  planar  secondary  face  to  be  twisted  is 
reduced,  and  the  L2  norms  for  sfl  and  sf2  become  identical.  For  both  sfl  and  sf2,  the  second  refinement 
(4x)  produces  a  reduction  in  the  L2  norm  of  about  3;  with  the  third  (8x),  it  increases  to  3.6.  The  reduction 
factor  for  sfO  is  about  5  for  all  refinements. 

The  y  flow  direction  results  for  test  2  arc  essentially  identical  for  all  three  shape  functions  (sfO,  sfl  and 
sf2);  this  is  because  the  secondary  faces  are  oriented  such  that  they  can  be  influenced  only  by  flow  in  the 
x  direction,  regardless  of  refinement.  To  see  this,  note  that  by  (44),  (56a),  and  (56b),  Wi*  =  Y2V(l/2)  x 
7j2x{  I  /2),  with  analogous  expressions  for  W2V  and  W2-;  the  Xt  =  Xi  —  Xo  vectors  are  differences  of  vectors 
that  have  the  same  x-component,  which  is  equal  to  the  distance  between  the  parallel  planes  normal  to  the 
A-axis,  hence  the  X2  vectors  lie  in  these  planes.  The  Y2  and  Z2  vectors  also  lie  in  these  planes;  hence,  all 
cross  products  in  the  W2  vectors  arc  parallel  to  the  A-axis.  Thus,  the  fluxes  fx 2  for  the  secondary  faces  arc 
null  for  for  either  sfl  or  sf2;  as  sfO  does  not  include  the  effects  of  a  secondary  flux  term,  all  three  results 
arc  identical.  For  all  three,  the  solution  degrades  with  the  first  refinement  (2x),  and  subsequent  refinements 
improve  the  L 2  norm  by  a  factor  of  3.3  to  3.7.  The  initial  degradation  in  the  L 2  norm  is  likely  caused  by 
non-planar  primary  cell  faces  in  the  refined  mesh. 

Test  3  examined  the  effect  of  non-planar  faces,  in  a  realistic  scenario  of  mesh  distortion,  on  the  L 2  norm; 
an  exact  solution  is  not  expected  under  these  circumstances.  The  2-D  mesh-face  distortion  can  be  measured 
by  the  sine  of  the  angle  between  poles  located  at  opposite  vertices  of  a  face;  the  average  sine  of  the  largest 
angle  between  poles  was  0.03  for  the  initial  mesh  (lx).  Table  I  indicates  that  the  flow  field  is  fairly  well 
resolved  for  the  initial  mesh,  regardless  of  the  shape  function  used;  sfl  and  sf2  are  only  marginally  better 
than  sfO.  Compared  to  test  1,  it  is  apparent  that  even  the  small  amount  of  facial  distortion  in  this  mesh 
degrades  the  solution,  masking  most  of  the  corrective  effect  of  the  quadratic  term  in  sfl  and  sf2.  With 
refinement  (8x),  the  average  sine  of  the  largest  angle  on  the  non-planar  surfaces  decreased  to  0.003;  the 
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Table  II.  Lr  norm  results  for  non-uniform  flow  tests.  Test  1:  homogeneous  medium;  comer-to-corner  flow.  Test  le: 
repeat  of  Lr  norm  calculation  for  test  1  with  exclusion  zone  around  flux  boundary  conditions.  Test  2:  heterogeneous 
medium,  reg:  orthogonal,  regular  mesh,  sf  1 :  shape  function  using  first  secondary  flux  correction  factor  in  conjunction 
with  the  quasi-random  mesh.  sf2:  shape  function  using  second  secondary  flux  correction  factor  in  conjunction  with 
the  quasi-random  mesh. 


test 

grid  refinement,  L2  norm 

no 

case 

lx 

2x 

4x 

i 

reg 

1.504  x  10“  3 

6.593  x  10“4 

3.436  x  10“4 

sfl,  sf2 

1.658  x  10“  3 

7.117  x  10“4 

3.704  x  10“4 

le 

reg 

1.109  x  10“  3 

2.407  x  10“4 

4.638  x  10“5 

sfl,  sf2 

1.114  x  10“  3 

2.478  x  10“4 

4.766  x  10“5 

2 

reg 

2.371  x  10“2 

1.424  x  10“2 

6.922  x  10“3 

sfl,  sf2 

2.489  x  10“  2 

1.473  x  10“2 

7.086  x  10“3 

decrease  was  uniform  and  first-order.  Concomitant  with  this  decrease,  the  reduction  factor  for  the  L 2  norm 
increased  from  2.6  to  3.6.  As  expected,  differences  between  results  from  sfO,  sfl  and  sf2  dissipate  as  the  role 
of  the  secondary  flux  dissipates  with  grid  refinement. 


9.  NON-UNIFORM  FLOW  TESTS 

In  this  section,  second-order  behavior  of  the  flux  truncation  error  (22)  is  examined,  as  well  as  non-uniform 
flow  in  heterogeneous  media.  For  all  tests,  a  regular  and  a  quasi-random  discretization  were  used;  for  the 
quasi-random  discretization,  all  x,  y  and  z  cell  faces  were  randomized  with  the  test  1  procedure  described  in 
the  previous  section.  Initial  grids,  before  refinement,  were  5x5x5;  the  domain  Q.  was  a  regular  cube.  For 
the  quasi-random  meshes,  non-linear  shape  functions  with  both  types  of  correction  factors,  sfl  and  sf2,  were 
tested. 

In  test  1,  with  homogeneous  K,  non-uniform  corner-to-corner  flow  was  obtained  by  injecting  a  unit 
flux  on  the  three  exterior  faces  of  cell  (1, 1, 1)  and  extracting  a  unit  flux  from  those  of  cell  (5,5,5).  With 
refinement  of  the  mesh,  the  flux  boundary  conditions  were  maintained,  preserving  the  original  facial  area 
of  the  applied  fluxes.  In  test  2,  a  heterogeneous  hydraulic  conductivity  field  was  generated  for  a  central  paid 
of  the  5  x  5  x  5  cell  mesh,  occupying  3x5x5  cells  of  the  mesh.  With  a  mean  value  of  about  100,  K  ranged 
over  three  orders  of  magnitude  (although  a  range  of  two  orders  was  most  common).  The  1x5x5  groups 
of  cells  at  each  end  of  the  heterogeneous  zone  were  assigned  a  uniform  K  of  100;  uniform  flux  boundary 
conditions  were  applied  over  the  exterior  faces  of  these  cells  to  create  a  mean  flow  in  the  ^-direction.  Grid 
refinement  preserved  the  structure  of  K.  In  these  simulations,  the  finest  (40  x  40  x  40)  discretization  was 
assumed  to  represent  the  “true”  solution;  L2  norms,  as  defined  in  (59),  were  obtained  by  comparing  coarser 
discretizations  with  the  finest.  Intermediate  discretizations  were  lOx  lOx  10  and  20x  20x  20.  The  L 2  norm 
results  for  these  simulations  are  summarized  in  Table  II. 

Note  that  the  errors  arc  generally  greater  than  those  of  Table  I;  this  likely  reflects  the  non-uniformity  of 
the  flow  rather  than  shape-function  issues.  That  Table  II  shows  the  same  results  for  sfl  and  sf2  indicates  that 
any  difference  between  them  is  overshadowed  by  the  errors  due  to  non-uniform  flow.  It  is  encouraging  that 
the  irregular-mesh  L2  norms  are  about  the  same  as  for  the  regular  mesh,  especially  for  the  unrefined  case, 
though  it  is  not  clear  which  aspects  of  the  treatment  of  mesh  distortion  are  responsible  for  this. 


final_report.tex;  6/11/2001;  9:04;  p.16 


17 


For  non-uniform  flow  in  a  homogeneous  medium  (test  1),  the  L2  norm  decreases  with  each  mesh  re¬ 
finement  by  a  factor  of  about  two,  indicating  first-order  convergence.  Though  (22)  is  second-order,  [1] 
demonstrates  that  convergence  will  be  degraded  if  there  are  singularities  in  the  exact  solution.  Flux  boundary 
conditions  on  parts  of  the  faces  of  c)Q  cause  singularities  at  the  points  where  the  nonzero  fluxes  meet  no-flow 
conditions.  In  test  le,  the  L?  norm  estimate  (59)  was  modified  to  exclude  errors  from  a  zone  surrounding 
each  of  the  flux  boundary  conditions;  then  the  L 2  norm  decreases  with  mesh  refinement  by  factors  varying 
from  4.5  to  5.2,  in  line  with  the  second-order  convergence  expected  from  (22).  For  the  heterogeneous  case 
(test  2),  the  L2  norm  results  produce  rates  of  improvement  which  vary  from  1.7  to  2.5  with  mesh  refinement, 
suggesting  first-order  convergence.  For  a  heterogeneous  problem  with  many  internal  discontinuities  in  K, 
this  is  expected,  as  these  discontinuities  can  result  in  singularities  in  the  exact  solution  [1], 


10.  CONCLUDING  REMARKS 

The  most  general  conclusion  of  this  investigation  is  that,  for  distorted  logically  rectangular  meshes,  the 
lowest-order  velocity  shape  functions  should  not  vary  linearly  with  space,  in  terms  of  either  the  velocity  or 
the  flux.  This  results  from  a  study  of  flux  conservation  in  a  general,  hexahedral  cell,  where  a  quadratic  term 
is  needed  in  order  to  match  the  flux  of  uniform  flow;  this  is  demonstrated  with  (18),  and  the  second-order 
estimate  of  the  error  of  quadratic  flux  interpolation  is  given  by  (22).  In  particular,  the  linearly  varying  flux 
of  the  Piola  transformation  cannot  match  uniform  flow  on  general  hexahedral  cells. 

This  quadratic  term,  however,  leads  to  an  additional  unknown  secondary  flux,  fx2,  which  must  be  in¬ 
dependently  estimated.  The  method  in  (5 1)— (52)  for  this  estimate  is  exact  for  uniform  flow;  in  general,  it 
is  at  least  concordant  with  the  second-order  flux  approximation  (18),  as  demonstrated  in  Appendix  D.  In 
addition,  we  have  found  that  weightings  for  (42)  such  as  (41)  and  (43),  which  adjust  the  velocity  magnitude 
by  comparing  a  differential  element  in  the  interior  with  a  differential  element  on  the  cell  face,  work  well  to 
interpolate  the  velocity  magnitude  to  the  interior. 

Meshes  with  non-planar  primary  and  secondary  cell  faces  will  lead  to  some  degradation  in  accuracy. 
Uniform  flow  tests  suggest  that,  as  long  as  this  type  of  cell  distortion  is  kept  minimal,  overall  accuracy 
is  not  greatly  affected.  In  particular,  the  L2  flux  error  for  a  mesh  containing  rather  smooth  non-planar 
secondary  faces  was  not  excessive.  Indeed,  if  at  least  two  faces  of  every  cell  can  be  made  parallel,  then 
shape  functions  which  compensate  for  twisted  planar-  faces  can  apparently  eliminate  this  error.  The  cost  is 
the  loss  of  generality  in  instituting  an  irregular-  grid.  Second-order  convergence  was  strongly  indicated  by 
the  reduction  factors  for  the  l 2  norms  for  all  the  uniform  flow  tests. 

In  regions  away  from  points  of  discontinuity,  the  non-uniform  flow  simulations  in  a  homogeneous  medi¬ 
um  also  show  second-order  convergence.  First-order  convergence  is  seen  in  non-uniform  flow  simulations 
in  heterogeneous  media.  Overall  convergence  of  most  practical  simulations,  which  will  have  singularities  or 
points  of  discontinuity,  should  be  first-order. 

The  results  reported  here  do  not  include  non-smooth  “random”  refinements,  in  which  angular  deviations 
do  not  decrease  (see  discussion  after  (23)).  Preliminary  computations  show  the  expected  first-order  conver¬ 
gence  in  many  cases,  and  less  in  some  others,  possibly  due  to  severe  distortions  caused  by  the  constrained 
procedure  that  generated  the  refined  meshes.  More  thorough  investigations  of  random  refinements  will  be 
reported  in  the  future. 


11.  APPENDIX  A 

A  form,  explicit  in  x,  for  the  logical  x  direction  cross  product  Y  x  Z,  is  derived;  the  logical  y  and  z  direction 
cross  products  are  similar.  From  (4)  and  definitions  of  the  covariant  vectors, 

Y  =  (1— x)Yo,  +  xYLv,  Z  =  (l-x)Z0x+xZlx,  (60) 
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where  Yo.v  =  Y|f=o,  Zov  =  Z|?=o,  Y \x  —  Y|x=i  and  Z\x  =  Zjx=i .  Manipulation  of  the  cross  product  gives  the 
desired  result: 

Y  x  Z  =  (1  —  x)2(Yox  x  Z0,)  +  x2(Yix  x  Zlx)  +  x(l  —  x)  [Y \x  x  Zox  + Yox  x  Zix] 

=  (1  —  x)(Y0v  x  Z0x)  +  x(Ylx  x  Zlx) -x(l-x)  [(Y0x-  YLv)  x  (Z0x-Zix)] .  (61) 


12.  APPENDIX  B 


By  (18)  and  (19)  (scaled  to  an  interval  of  length  Ax,  with  x/Ax  in  the  role  of  x),  along  with  (21), 

fx(x)  -  Qfx{x)  =  G(x, ^)/"(4)  d%  +  y  ( i  -  y  )  fix  ^y  ^  .  (62) 

As  the  area  bounded  by  the  graph  of  G(x,  •)  is  a  triangle  of  base  Ax  and  (negative)  height  x(x/Ax  —  1), 


L ' ■ 4  =  'ixAx  (£  -  0 = (‘  -  £)  ■ 

Combining  (62)  with  (63),  we  obtain 


r  Ax 

fx(x)  -  Qfx(x)  =  /  G(x. 

Jo 


(ft,. 


(63) 


(64) 


We  evaluate  the  bracketed  term  in  (64),  differentiating  fx  twice,  where  the  appropriate  analogue  of  (16)  is 

AW  =  ('  -  s) A  w ~  b  (‘  -  £) (65) 


This  gives 


m  -  (y)  =  + (i  -|)  A(B + ^ 


flxi^,)  —  fix  (  y 


+s  (s  -  s)  A®  - 1  (■  ■ -  h) + s^® + s*®  ■ 


(66) 


The  expressions  in  (17a)-(17c)  and  some  rearrangement  then  yield  (22). 


13.  APPENDIX  C 

To  verify  the  first  relation  in  (26),  define  the  subelement  Qxc  of  Q,  consisting  of  points  r(x,y,z)  for  0  <  x  <  c, 
0  <y  <  1,  0  <  £  <  1,  and  the  linear  mapping  Lxc(x,y,z)  =  (cx.y.z).  Then  roLxc  maps  Q  to  Qxc,  roLxc  is 
trilinear,  and  Fxc  is  on  the  boundary  (i.e.,  is  a  primary  face)  of  Qxc.  Taking  vxc (x. y, z)  =  (cx,  0,0)  on  Q,  we 
apply  the  Piola  transformation  (24)  to  r  o  Lxc  and  r  and  obtain 

vxc(roLxc(r))  =  (7(roLxc))(f)_1D(roLxc)(f)vxc(r)  =  (/(roLxc))(r)_1  D(roLxc)(f)(cx,0,0)r 

=  [(T(r))(Lxc(r))c]_1Dr(Lxc(f))(c2x,0,0)7'  =  (/(r))(Lxc(f))~'  Dr(Lxc(r))(cx,0,0)7’ 

=  (7(r))(Lxc(r))^‘  Dr(Lxc(r))v(Lxc(r))  =  v(r(Lxc(r)))  =  v(roLxc(r)).  (67) 

Thus,  the  transformed  vector-valued  functions  vxc  and  v  take  the  same  values  on  Qxc,  and  hence  on  Fxc.  Let 
s  denote  the  variable  on  Fxc,  6  on  Fx\ ,  and  s  on  Fxc,  with  s  =  r(s) ,  s  =  Lxc (d),j  =  ro  Lxc (a) .  Now,  because 
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Fxc  is  on  the  boundary  of  Qc,  we  can  use  (67),  (25),  vxc  ■  n  =  c  on  Fx\,  and  the  fact  that  the  surface  Jacobian 
of  Lxc  is  1  on  Fx\  to  obtain 

[  \-n{s)z{s)ds  =  [  \xc  -n(s)z(s)ds  =  [  \xc-h{a)z{roLxc(a))da 
7  Fxc  7  Fxc  7  Fx\ 

=  I  cz(roLxc(o))dd  =  f  cz{r{s))ds ,  (68) 

7  Fxl  7  Fxc 

as  desired. 

For  the  second  relation  in  (26),  similarly  take  Qyc  to  be  those  r(x.y.z)  for  0  <  jt  <  1,  0  <  y  <  c,  0  <  £  <  1, 
L yc(x,y,£)  =  (x,cy,£),  and  vyc(x,y,£)  =  (x,0,0)  on  Q.  Then 

cvw(roLvc(f))  =  (/(r  o  Lyc) )  (r)  ~ 1  cD  (r  o  Lyc)  (r)  vw  (r )  =  (/(roLvc))(f)_1  cD(roL>IC)(f)(f,0,0)7’ 

=  [(7(r))(L>r(f))c]  1  cDr(Lvc(r) )  (x,  0, 0) T  =  (J(r))(Lyc(r))  1  Dr(Lvc(r))(x:0,0)T 
=  (/(r))(Lyc(f))^1Dr(L3,c(f))v(L3;c(f))  =  v(r(Lvc(r)))  =  v(roL>,c(f)).  (69) 

An  ai'gument  like  (68)  then  yields  the  second  relation  in  (26). 


14.  APPENDIX  D 

We  show  here  that  the  secondary  flux  approximation  (51),  which  completes  the  definition  of  our  velocity 
shape  functions,  does  not  affect  the  order  of  accuracy  of  quadratic  flux  interpolation.  From  (17c),  (18),  and 
(46),  the  exact  secondary  flux  to  be  approximated  by  fx 2  is 

fxl  =  At (1/2)  =  [ 1  [ 1  q(l/2,y,£)  •  (Y2,(£)  x  Z 2x(y))dydt  (70) 

Jo  Jo 

As  shown  in  (64),  scaled  to  x-dimension  Ax,  the  contribution  of  fxi  to  the  error  (fx  —  Qfx )  (x)  in  the  quadratic 
flux  interpolation  in  (22)  is 


rAx 

LfAx)  -  Qfx{x)  =  /  G(x,^) 

Jo 


Ax1-1*2 


df. 


(71) 


In  (71),  integration  produces  a  factor  of  Ax  and  G(x, q)  =  O(Ax);  these  factors  are  canceled  by  —2/Ax2. 
Thus,  to  preserve  the  second-order  accuracy  of  (22),  we  seek 


fxl- fxl  =  0(Ax2). 


(72) 


If  (72)  holds,  then  the  interpolated  fluxes  Qfx  associated  with  the  velocity  approximation  Vr  of  q  in  (46)-(50) 
will  maintain  second-order  accuracy. 

From  (51),  (52),  and  (70),  we  see  that 


fxl  -  fxl 


■  (Y2x  x  Z2,)  dydl 


(73) 


Assuming  a  smooth  refinement  as  in  the  discussion  following  (23),  Y2_t  and  Z2v  are  0( Ax),  hence  Y2a  x  Z2a-  = 
0(Ax2).  The  sum  in  (73)  involves  discretization  in  all  three  directions,  so  that  we  must  consider  refinement 
in  y  and  £  as  well  as  x.  We  seek  the  error  on  the  x-face,  which  is  refined  in  reference  space  into  1/AyA z 
subfaces  of  v-dimension  Av  and  £-dimension  Az.  By  keeping  the  same  relationship  between  reference  and 
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coarse-cell  coordinates,  independent  of  refinement,  we  keep  the  same  X,  Y,  and  Z  vectors  and  the  same 
Jacobian  J  (otherwise,  we  would  have  to  rescale).  Thus,  on  one  such  subface,  we  seek 


f, 


x2  r 


nAy  J  1  /  fa  \ 

-  Yj  Yc(iAxJAy,kAz)  -q[—,y,z) 

8  feo  V2  J 


L  774=  0 
=  0(Ax2AyAz)1 


■  (Y2,  x  Z2x)  dydz 


(74) 


so  that  the  sum  of  the  errors  on  the  subfaces  will  remain  0(Ax2).  The  estimate  that  we  prove  below  is 


fx2r  -  fxlr  =  0(( Ax2  +  Ay2  +  Az2)  Ar2AyAz) , 


(75) 


which  is  a  much  stronger  result,  showing  that  the  error  of  the  secondary  flux  approximation  is  negligible 
compared  to  that  of  quadratic  flux  interpolation. 

Split  (74)  into  two  terms  7)  and  7),  replacing  Vc  with  q  in  7).  Apply  the  midpoint  integration  rule  in  y 
and  z  in  T\  to  obtain 


Ti  = 


nAy 

1 


l  Y  q('Ar,  j  Ay,  kAz)  -q(^-,y,z 
L°  774=0  \Z 


(Y2jc  X  Z2t)  dydz 


Y2x[^)xZ1x(1- 


AyAz 


+  0((Ay2  +  Az2)Ar2AyAz)  =  0((Ar2  +  Ay2  +  Az2)Ar2AyAz), 


(76) 


where  the  last  equality  estimates  the  difference  between  the  sum  and  the  midpoint  value  of  q. 
For  the  remaining  term  of  (74),  we  have 


T2  = 


Az  r  A  y 


n 

Jo  Jo 


1 


Y  (Vc-q)(/Ax,;Ay,kAz) 
L  U  774=  0 

=  {Eox  +  E\x  +  Eoy  +  Ely  +  E0z  +  E\z)  /2, 
where  the  E  terms  arc  associated  with  the  cell  faces;  for  example, 


(Yzr  x  Z2r)  dydz 


Eox  — 


A z  rAy 


rixz  r 

Jo  Jo 


T  Y  (' Vx  -  qx)  (0,  ./Ay,  kAz) 

L  74=0 


(Y^xZ  2x)dydz. 


(77) 


(78) 


In  (78),  qx  and  Yx  are  multiples  of  X,  as  in  (32)  and  (50).  By  construction,  with  property  (UE2)  and  (46), 

Vx  -  -jfi.xofx o  -  7  Av(0) 


is  defined  such  that  the  normal  component 


„  _  1X-(YxZ)|YxZ|,  fM 

r  \-wr..rrl  a  /^\^  JX 0 


J  |Y  x  Z|  Ar(0)  4,(0) 

is  constant  on  the  face.  By  (79),  Yx  satisfies 

Y  (yx  ■  ni)  (0,  ./Ay,  kAz )  = 


(79) 


74=0 


4,(0) 


(80) 
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Now  by  (14)  and  the  midpoint  integration  rule 

fAz  Mv 


ray  ~  » 

/vO  =  /  /  qx-(YxZ)dyd£  =  (qx  •  (Y  x  Z))(0,  Ay/2, Az/2)AyAz+ 0(AyAz(Ay2  +  Az2)), 

J  0  Jo 

Ax(0)  =  r  r  |Y  x  Z| dydz  =  |Y  x  Z|(0, Ay/2, Az/2)AyAz  +  0(AyAz(Ay2  +  Az2)), 

Jo  Jo 


fx  0 

A,(0) 


(qx-  (YxZ))(0, Ay/2, Az/2) 


|YxZ|(0,  Ay/2,  Az/2) 


+  0(Ay2  +  Az2)  =  (qx  •  nf )  (0,  Ay/2,  Az/2)  +  0(Ay2  +  Az2) 


1 


t 


=  t  Y  (<1x  '  n-?)  (°’  AAz)  +  0(Ay2  +  Az2) . 
4 14=0 


By  (80)  and  (81), 


Now 


T  I((Vx-qx)-n,) (0, jAy, ^Az)  =  0(Ay2  +  Az2) . 


4 14=0 


/•Az  /•A); 

/  /  (Vx  —  qx)  •  n.f  |Y  x  Z|(0, y,z)  dydz  =  fxo  —  fxo  =  0, 

Jo  Jo 

hence  for  some  (yo,£o)  we  have 

0=  ((Vx-qx)-n.v)(y0,£o)  =  (|VX  -  qx|  j^|  •  n.f)(y0,Zo), 

whence  (Vx  —  qx)(vo,£o)  =  0,  which  implies  that 

|Vx-qx|  =  0(Ay  +  Az). 

Because  n_v(0,yAy,Mz)  =  n^(0,  Ay/2,  Az/2)  +  0(Ay  + Az),  (82)  and  (83)  give 

t  £  (Vx-qx)(0,  ./Ay,  Mz)  •  nf  (0,  Ay/2,  Az/2) 

4 14=o 

=  0(Ay2  +  Az2)  +  0(Ay  +  Az)  0{  Ay  +  Az)  =  0(Ay2  +  Az2) , 


(81) 


(82) 


(83) 


from  which  we  immediately  obtain 

t  L  |Vx-qx|(0,iAy,/cAz)  =  0(Ay2+Az2). 

4 14=o 


Then  (78)  and  (84)  yield 

Eqx  =  0((Ay2  +  Az2)Ar2AyAz), 
and  a  similar  argument  for  the  other  five  terms  in  (77)  leads  to 

T2  =  0{  (Ax2  +  Ay2  +  Az2)  Ar2AyAz) . 


(84) 


(85) 


Finally,  from  (76)  and  (85),  we  obtain  (75),  whence 

fx 2  -  fx 2  =  0((  Ax2  +  Ay2  +  Az2)Ax2),  (86) 

as  desired.  This  is,  in  fact,  a  significantly  stronger  result  than  (72),  showing  that  the  error  in  the  secondary 
flux  approximation  is  negligible. 
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If  the  refinement  is  not  smooth  (“random”  refinement),  then  we  find  that 

fx 2  -  fx2  =  0(Ar  +  Ay  +  A z) ,  (87) 

so  that  the  first-order  error  in  the  discussion  after  (23)  is  preserved.  To  see  this,  trace  through  the  argument 
above  and  observe  what  is  lost.  The  factor  Ax2  is  lost  throughout  because  Y^  x  Z^r  =  0(1)  instead  of 
0(Ax2).  One  order  is  lost  in  every  estimate  involving  the  midpoint  rule,  because  Y  x  Z  or  n;  can  have 
derivatives  of  0(Ay_l  +  Az~!).  One  order  is  lost  in  (84),  because  (83)  still  holds,  but  rif  may  vary  by  0(1). 
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